Diffusion of curcumin in PLGA-based carriers for drug delivery: a molecular dynamics study

Context: The rapid growth and diversification of drug delivery systems have been significantly supported by advancements in micro- and nano-technologies, alongside the adoption of biodegradable polymeric materials like poly(lactic-co-glycolic acid) (PLGA) as microcarriers. These developments aim to reduce toxicity and enhance target specificity in drug delivery. The use of in silico methods, particularly molecular dynamics (MD) simulations, has emerged as a pivotal tool for predicting the dynamics of species within these systems. This approach aids in investigating drug delivery mechanisms, thereby reducing the costs associated with design and prototyping. In this study, we focus on elucidating the diffusion mechanisms in curcumin-loaded PLGA particles, which are critical for optimizing drug release and efficacy in therapeutic applications. Methods: We utilized MD to explore the diffusion behavior of curcumin in PLGA drug delivery systems. The simulations, executed with GROMACS, modeled curcumin molecules in a representative volume element of PLGA chains and water, referencing molecular structures from the Protein Data Bank and employing the CHARMM force field. We generated PLGA chains of varying lengths using the Polymer Modeler tool and arranged them in a bulk-like environment with Packmol. The simulation protocol included steps for energy minimization, T and p equilibration, and calculation of the isotropic diffusion coefficient from the mean square displacement. The Taguchi method was applied to assess the effects of hydration level, PLGA chain length, and density on diffusion. Results: Our results provide insight into the influence of PLGA chain length, hydration level, and polymer density on the diffusion coefficient of curcumin, offering a mechanistic understanding for the design of efficient drug delivery systems. The sensitivity analysis obtained through the Taguchi method identified hydration level and PLGA density as the most significant input parameters affecting curcumin diffusion, while the effect of PLGA chain length was negligible within the simulated range. We provided a regression equation capable to accurately fit MD results. The regression equation suggests that increases in hydration level and PLGA density result in a decrease in the diffusion coefficient. Supplementary Information The online version contains supplementary material available at 10.1007/s00894-024-06023-x.

consuming simulations such as MD ones.Furthermore, the expected experiments result to be often redundant, adding little or even no new information.Hence, TM results particularly adapt to this case, since it does not require too many trials.An orthogonal array (OA) is used to reduce the number of simulations and obtain reasonable information. 2Writing the orthogonal array matrix involves a procedure also known as the fractional factorial design technique.
In the considered domain, the three input variables to be explored are: the hydration level, the PLGA length, and density.Considering these three variables and three levels per each one of them, the TM required a number of trials equal to 3 2 .Hence a L9 matrix was designed, meaning that the design space is composed by three values for each sample and thus nine possible configurations to be simulated.
After the DOE structure and the orthogonal matrix have been written, the next step is to conduct the matrix experiment.The following protocol was followed to build the target MD configurations: 1. Placing the curcumin in the center of a sphere formed by packed PLGA chains with a given length.
2. Loading the packed ensemble in an opportunely large cubic box.
3. Performing an energy minimization of the system to relax the compacted chains.
4. Introducing a number of water molecules in the computational box needed to reach the target hydration level, performing a further energy minimization after the solvation process.The water density was computed with respect to the solvent accessible volume in the box. 5. Performing a series of NVT and NPT simulations to equilibrate the temperature and pressure of the system and reach the target PLGA density, while adjusting the number of water molecules with the updated solvent accessible volume in the box.An accept-able margin of error on the final PLGA density was established to be ±5% with respect to the target values.
The final box size, number of water molecules and number of total atoms for DOE1 and DOE2 are summarized in Tabs.S3 and S4.After performing this preparation procedure for all configurations prescribed by TM, NVT production runs are carried out to generate statistically relevant trajectories (5 ns) and finally to post-process the diffusivities of water, PLGA, and curcumin molecules.

Supplementary Note 3: Simulation convergence
The time evolution of the MSD for all the simulations of DOE1 and DOE2 are presented in Figs.S1, S2, S3, S4, S5 and S6.Each D value was calculated from the linear fitting of these MSD trajectories, following the Einstein relation.As it can be noticed from Figs. S1, S2, S3, S4, S5 and S6, all simulations are stable, with diffusivity values that pass the Tukey's method check, proving the convergence of results.
Additionally, a test simulation (i.e., sim 7 of DOE1) was conducted for up to 10 ns to check the system's convergence.As shown in Fig. S7, the mean squared displacement (MSD) of each species remains stable over time.We calculated the diffusion coefficients for this extended simulation, obtaining the following results: for water, D = 11.176• 10 −12 m 2 /s; for curcumin, D = 0.649 • 10 −12 m 2 /s; and for PLGA, D = 0.592 • 10 −12 m 2 /s.The percentage relative errors between the 10 ns and 5 ns simulations were then computed for each species, yielding a discrepancy of 4.82% for water, 2.14% for curcumin, and 0.42% for PLGA.Thus, we considered the differences between the diffusion coefficients to be negligible and regarded all the 5 ns simulations as sufficiently long to guarantee convergence.The resulting diffusion coefficients and the linear fitting interval times for each simulation are summarized in Tabs.S1 and S2 for DOE1 and DOE2, respectively.

Figure S1 :Figure S2 :Figure S3 :Figure S4 :Figure S5 :Figure S6 : 2 )
Figure S1: Time evolution of solvent MSD in DOE1 with linear fitting for D calculation, from configuration 1 to 9 (from left to right; from top to bottom).

Table S2 :
Diffusion coefficients and linear fitting interval time for DOE2.N.A. refers to Not Applicable, when water molecules are not included in the simulation.−12 m 2 /s] Linear fitting time interval [ps] D [10 −12 m 2 /s] Linear fitting time interval [ps] D [10 −12 m 2 /s] Linear fitting time interval [ps]

Table S3 :
Box size, number of water molecules and total number of atoms for DOE1.

Table S4 :
Box size, number of water molecules and total number of atoms for DOE2.